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Abstract 

The extinction transition in the presence of a localized quenched defect is studied numerically. 
When the bulk is at criticality, the correlation length diverges and even an infinite system cannot 
"decouple" from the defect. The results presented here suggest that, in 1+1 dimensions, the 
critical exponent 5 that controls the asymptotic power-law decay depends on the strength of the 
local perturbation. On the other hand, the exponent was found to be independent of the local 
arrangement of the defect. In higher dimensions the defect seems to induce a transient behavior 
that decays algebraically in time. 
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I. INTRODUCTION 



In a recent article, the results of a spatial-genetic study of honeybees in the Sahara 
desert []| have been published. The authors concluded that these bees (Apis mellifera) 
first arrived to the desert around 10,000 years ago, when the Sahara was a green savannah. 
Although the contemporary Sahara is inhospitable to these bees, they do survive in many 
oases that litter the desert. The same scenario holds for many existing communities of 
animals: their habitat is made of relatively small patches, and the sustainability of the whole 
system is determined by the rates of local extinction and recolonization. When migration 
among "oases" is relatively weak, such a system is known as a metapopulation j^j]. 

In the study of metapopulation persistence, many researchers are using the stochastic 
patch occupancy models (SPOM) 3|-[5| as a paradigmatic framework. These models assume 
that each habitat patch is either occupied or empty, i.e., that the timescale in which the 
population reaches its local carrying capacity equilibrium is relatively small. Neglecting the 
details of the in-site dynamics the theory is then focused on extinction and recolonization 
events. In the language of stochastic dynamics, we are dealing with a contact process 
taking place in heterogenous environment. While the properties of the contact process and 
its generalizations are quite well understood when it takes place on a uniform lattice, much 
less is known in the presence of a quenched disorder. In this paper we are dealing with a very 
simple (generalized) contact process perturbed by a local heterogeneity. Interestingly 

it turns out that even a local defect may alter the behavior of the whole system at the most 
important parameter regime, close to the extinction transition. 

On a homogenous lattice the contact process is known to admit an extinction transition 
when the birth rate is below some (nonuniversal) threshold. For any extinction transition 
on a homogeneous substrate and a single absorbing state, Grassberger jgj] and Janssen {^J 
conjectured that the microscopic details of the stochastic process are irrelevant close to the 
extinction point and the transition belongs to the directed percolation (DP) universality 
class. The basic rationale behind this conjecture is that a spatially extended system decom- 
poses, close to the transition, into active and inactive zones, where after each typical period 
of time there is certain probability for an active state to die, to survive, or to infect its inac- 
tive neighbors. If these regions are considered as lattice points onarf dimensional array, the 
chance of an active site to survive or to infect its neighbors within a unit time is equivalent 
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to the chance that a bond exists between a lattice point at time t and its neighbors in a 
subsequent replica of the system at t + 1. Accordingly, the extinction transition happens 
when the bond density is exactly at threshold for an infinite cluster in a d + 1 dimensional 
system, and the transition belongs to the directed percolation universality class in d + 1 
dimensions. 

The Grassbrger-Janssen conjecture has proven to be extremely robust, and a large number 
of stochastic models that admit an extinction transition were shown to belong to the DP 
equivalence class if the substrate is homogenous 
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11] . It was further shown that spatio- 



temporal substrate noise (i.e., birth-death rates that fluctuate in space and time with only 
short range correlations) is an irrelevant perturbation close to the transition, so small noise 



is averaged out and leaves the DP transition unaffected [10J, although it can change the 
location of the transition point [13L The case of diffusively correlated disorder 
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17] was 



considered recently by Dickman [18| and by Evron et. al. [19] with somewhat different 
results: while [isl reported that the critical exponents differ from those of DP transition, in 
191 ] it was suggested that only the off-transition behavior of the system is governed by local 
adaptation to favored regions, since the scaling function describing the behavior away from 
the transition shows significant deviations from the known DP behavior. 

In practice, however, it is hard to avoid quenched heterogeneities from the substrate on 
which the process takes place, and this in true in particular for metapopulations [2J or the 



spread of an epidemic 
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211 ] . It turns out that quenched (time independent) disorder is a 



relevant perturbation 12[ and seems to change the nature of the transition. 

When the system is subject to a global disorder (the spatial range in which the hetero- 
geneity occurs is not compact) a Griffiths phase exists between the active and the inactive 



regions 



141 ] . In the parameter region that corresponds to the Griffiths phase the survival 



of an occupied patch depends on the local properties of the substrate, not on activation by 
neighboring regions. In particular, for each time scale the process stays alive due to the 
spatial domains that admit high carrying capacity [yj, lj|. Although stochastic fluctuations 



guarantee extinction for any localized active patch in the absence of migration 22M24J] . the 



time to extinction grows exponentially with the carrying capacity of the spatial domain. 
This implies that exponentially rare spatial patches, with high birth rate, support the pop- 
ulation for exponentially long times. An optimization argument [l^] shows that in such a 
case the survival until t is dominated by rare spatial fluctuations of linear size L ~ log(t); 
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accordingly, the density falls algebraically with time. The Griffiths phase is located between 
an extinction region, where essentially no oases exist, and the active phase, where good oases 
infect each other to yield a never dying process. 

In this paper we intend to focus on a much simpler problem: the contact process in the 
presence of a localized (compact) spatial heterogeneity, like an oasis in the desert or an 
ecological trap (sink habitat). When the bulk system is off criticality it is segregated into 
effective patches of size £, and the effect of a compact inhomogeneity extends only to this 
distance, so it cannot change the long-time properties of an infinite sample. The situation 
changes when the bulk system is at criticality: in that case the correlation length decay 
obeys a power-law, and in principle a single oasis (or sink) may affect the properties of the 
bulk. Here we present numerical evidences suggesting that this is indeed the case, at least 
in one dimension. A similar phenomenon characterizes the physics of the zero temperature 
Kondo effect, when a single magnetic impurity changes the behavior of an infinite system as 
it dictated the density of states close to the Fermi level where correlation length diverges. 

Before proceeding to the results, let us describe the numerical procedure. We have 
simulated a spatial version of the SIS model for epidemics, which is a simple generalization 
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On each site 



of the contact process, utilizing the simulation technique presented in 
there are N agents, from which N — I are susceptible and / are infected. If the infection 
rate is a and the recovery rate is /3, the chance of an infected agent to recover within At 
is (3At. On a one dimensional lattice a susceptible agent on the n-th site has I n infected 
neighbors in its own site and I n ±i in the neighboring sites. If the chance to infect a local 
susceptible is a(l — x)/N and the chance to infect a neighbor susceptible is ax/(2iV), 
the average number of infections on the n-th site within At is given by q\(N — I n ) where 
qi = Ata [(1 — x)/„ + X-fn+i/2 + X-^n-i/2] and the average number of recoveries is q2l n where 
q 2 = Atf3I n . Once these numbers were calculated for all the chain, the system is updated 
by drawing, for any site, two numbers from a binomial distribution, one with an average 
qi(N — I — n) and the other with an average q2l n - In{t) is than updated to I n (t + At) = 
I n {t) + Biql, N — I n ) + B{q2, I n ). This procedure may be considered as a stochastic Euler 
integration of the system dynamics, and it converges to the pure contact process in the limit 
N = 1, x — 1/2- Without loss of generality we have taken = 1 in all cases considered 
below, and the generalization to higher dimensions is trivial. 



4 



II. DEFECT DEPENDENT DECAY EXPONENT IN ONE DIMENSION 



Let us start presenting our results for the contact process, N = 1. The critical birth 
(infection) rate in that case is known to be a c = 3.2978. For a homogenous system at 
criticality the survival probability, i.e., the chance of a cluster, grown from a single seed, to 
be active after t steps, is given by P(t) = t~ s , where the Id value is So ~ 0.159. 

Figured] shows what happens when a single site for which a = a c +5a has been introduced 
into a system which is otherwise at criticality (see schematic sketch in the upper-right part of 
the figure). When the process is ignited with a single particle at the heterogeneity (arrow), 
P(t) is still a straight line on a double logarithmic plot, but the slope depends on Aa: for 
source (oasis) the slope is smaller than 5 an d for sink the decay is faster. This indicate 
the failure of the critical system to decouple from the defect when the correlation length 
diverges. Figure |2] shows 5(Aa)/5 between a perfect sink and a very strong source. 

What happens when the process is ignited away from the heterogeneity? If the oasis/sink 
is located at the origin and the process starts when a seed (a single infected person) is 
positioned at x, P(t) must decay with the exponent 5 until the system starts to "feel" the 
oasis, then cross over to the appropriate value of 5. When the system is at criticality the 
spread of a perturbation is controlled by the ratio between the spatial and the temporal 
correlation length, i.e., a perturbation is traveling a distance x within a time t ~ x z , where 
for a Id DP transition z = 1.58. This implies that after a time t the source/sink affects a 
distance x ~ t om . 

Figure E] illustrates the effect of the distance from the transition point. Here the seed 
was deposited at Ax = 15 where the origin is an oasis, a sink or a regular site for which 
a = a c . P(t) for all cases is identical until the system start to get some information about 
the structure at the origin, than it splits into three different slopes (left panel). In the right 
panel there is a sink at the origin and P(t) is shown for different distances of the seed from 
the heterogeneity: one sees that the short time decay has the exponent of a homogenous 
system, the long time decay also have the same slope, but the crossover time depends on 
the distance. 

The same behavior is depicted in Fig. |H Here the local particle density p(x,t), starting 
from a fully occupied lattice, is shown: a snapshot at t = 1000 for different types of special 
site (inset) and consecutive time frames for a system with a sink. Since at the close region 
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FIG. 1: The survival probability P(t) vs. t on a double logarithmic plot for a system of L=5000 
sites with periodic boundary conditions. Only a single particle is allowed on each site, so the 
process is equivalent to the contact process. At the bulk the infection rate is a c = 3.2978, and only 
on the "exceptional" site at the origin a = a c + Act. When the system is homogenous the slope 
is So ~ 0.16, as expected from the theory. If on the central site a < a c (sink) the slope is larger, 
and for an oasis (a > a c ) it is smaller. The results here reflect an average over 10 6 processes, all 
ignited with a single agent on the source/sink site. 



p(t) decays like t s ( Aa ) and in the far region pit) ~ t So , on a double logarithmic scale a 
straight line appears for p{x) in the intermediate region, with a slope 

logjpot- 5 ^) - logjpot- 50 ) J(Aa) - 6 

s = 2 z • (1) 

This slope agrees with the measurements from our numerical experiments. 



III. UNIVERSAL RESPONSE TO LOCAL HETEROGENEITY 

Surprisingly, our numerics suggest that the effect of a localized heterogeneity is also 
universal, i.e., it is independent of the spatial structure of the source/sink and of N, the 
carrying capacity of a site. This feature is demonstrated in Figure |5j here the heterogeneity 
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FIG. 2: The decay exponent 5(Aa) vs. the strength of the heterogeneity Aa, from a sink where 
the birth rate is zero (Aa = — a c , i.e., a = 0) to a very strong oasis. 

is made of two exceptional sites, and the distance between the two, 2d, is a free parameter. 
When the process starts with a seed between the two sites, the short time decay of P(t) 
is, of course, nonuniversal. However, after the process "probes" the two sites the long time 
behavior is universal and independent of the distance between the sinks. We have carried 
out different numerical experiments, with two or three, equal or unequal, exceptional sites 
(and for two sources, source and sink, two sinks with different as etc.). In all of them the 
tail exponent is independent of the spatial structure of the compact heterogeneity. 

The same observation holds if one allows more than one agent at any site. In Fig. [6] 
P(t) is plotted for N = 50. There is a transient region in which the system followed the 
deterministic decay rate P ~ 1/t (as suggested by the equation / = a^V 2 / — al 2 , which is 



the deterministic limit of the stochastic process denned above 28[). However, beyond this 
transient region the stochastic exponent controls the decay, with different exponents for the 
homogenous and the inhomogeneous case. The values of the asymptotic decay 5 seems to 
be very close for N = 1, N = 10 (not shown) and A^ = 50 (differences smaller than 0.02, 
which is more or less the numerical error in the measurement of the slope), so it seems that 
the exponent is independent of A" either. 
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FIG. 3: Left panel: P(t) vs. t (loglog scale) for a contact process that starts with one agent at a 
distance Ax = 15 from the heterogenous site. Three cases are compared: an oasis (Aa = +7), a 
sink (Aa = —2.5) and a homogenous system. Right panel: same, with a sink (Aa = —2.5) at the 
origin, where the three lines correspond to different Ax. 



IV. TWO AND THREE SPATIAL DIMENSIONS: CONVERGENCE TO THE HO- 
MOGENEOUS CRITICAL BEHAVIOR 

All the results presented so far are for a Id system. In higher dimensions the situation 
is different: here the effect of a compact heterogenous source/sink manifests itself only as 
a transient, but the asymptotic behavior of P(t) follows the exponent So that characterizes 
the homogenous system, as shown in Figures [7J Careful examination of the results seems to 
indicate that the long time decay is a combination of two exponent of the form: 

Pit) ~r <5 ° (i + At~ Sl ) (2) 

for a sink and 

for a source, i.e., that the transient behavior itself decays algebraically in time, with an 
exponent larger than So- 

The qualitative difference between the behavior in Id and in higher dimension suggest 
that the effect of a single defect has to do with the chance of a random walker to return 
to the origin. At criticality the average spatial size of the contact process is not shrinking 
or growing, and since there is no preferred direction the process visits spatial points at 
random (this feature has been already considered in the theory of persistence time 
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FIG. 4: The main panel shows a few "snapshots" of the spatial density profile, p(x), for a contact 
process that starts from a fully occupied lattice with a sink (Aa = — a c ) at the origin. Each line 
reflects the average results of 1000 histories. While in the far region the density decays algebraically 
with 5q, in the close region the decay is faster. The region affected by the sink is also growing. The 
inset shows the density profile at fixed time, t = 1000, for different types of heterogeneity at the 
origin. Clearly, in the far region there is no difference, and the distance in which the sink/source 
affects the system is independent of its strength and character. 

Since in above 2d there is a good chance that the process never returns to the oasis, and 2d 
is marginal, only in Id the oasis changes the critical exponent. 

V. DISCUSSION 

Along this paper we have dealt with a single oasis/sink coupled to a critical bulk. This 
problem is quite similar to the one recently considered by Barato and Hinrichsen [26|, in 
which an inert spatial system (at the bulk there are no birth-death processes, only diffusion) 
is coupled to a single active site on the boundary. In this setup the active site admits only 
one agent (N = 1) with and is "critical" in the sense that the birth rate on this site is equal 
to the death rate. These authors reported a decay exponent 1/6 for that system. Later, 
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FIG. 5: The survival probability P(t) vs. t for a contact process. Here the bulk is at criticality 
with two, unequal, sink sites, one with a = a c — 1 and for the other a = a c — 2. The process starts 
when a single seed is introduced at the origin, which is the midpoint between the two sinks, with 
a distance d from the seed to any of the heterogeneities. As d increases, the time it takes to the 
process to interact with the sinks is growing, still the asymptotic behavior is the same (the slope 
of all these lines is the same within the measurement error). 



Burov and Kessler [27| have solved analytically the problem for N = oo and found 5 = 1/4; 
For any finite iV their numerics indicated algebraic decay with different exponent. 



In some sense, the model considered here and the model of [261, [27J are similar, as they 
both deal with an exceptional site coupled to a critical bulk. The main difference is that in 
our model the bulk is active and the criticality is "quantum", while the model considered 
before has an inert bulk. Still it will be interesting to check the response of the inert bulk 
model to a spatially structured defect. 

An interesting question that may be asked in general with regard to these systems is about 
their deterministic limit, i.e., the behavior when A" — > oo. For the deterministic system the 
critical point is at a = 1 and a single defect coupled to a critical medium satisfies an equation 
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FIG. 6: The density pit) vs. time for a system with carrying capacity of 50 agents per site. With 
X = 0.2, the transition for the homogenous system has been identified at a c = 1.1875. At criticality 
(black line) the decay starts with exponent 5 = 1, as expected from the mean-field equations, then 
when the stochasticity becomes important it crosses over to 5q of the DP transition (the dashed 
line is the linear fit for the tail). On the other hand, if a sink (site with a = 0) is introduced at 
the origin the long-time exponent is larger (red line) and the best fit to the tail (dashed) yields 
5 = 0.41. The results reflect an average over 10 6 histories. 



like 

I = V 2 / + a(x)I — J 2 , 

where a(x) is takes nonzero values at the spatial region associated with the defect. Since 
any potential well in Id supports a localized state, so does the linearized evolution operator 
£ = V 2 + a(x). Thus for any oasis, with any spatial structure, the system must admit 
a localized colony (its amplitude is dictated by the nonlinear term) that lives forever 29]. 



P(t, N) is thus singular at the deterministic limit. In higher dimensions, on the other hand, 
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FIG. 7: The survival probability for a process initiated with a single seed at the origin for 2d (a) 
and 3d (b). The site at the origin has different reproduction rate. Either for source or for sink, the 
long-time decay seems identical and converge to the known values for the DP exponents in these 
dimensions. The bulk values of a c are 1.6488 for 2d and 1.3168 for 3d. 

not any oasis supports a localized state, so the deterministic limit may be non-singular. 
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